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ABSTRACT 


An  experimental  and  theoretical  analysis  of  a  light  shallow 
clamped  arch  under  a  central  concentrated  load  is  performed.  Arches 
covering  a  wide  range  of  arch  geometry  are  considered.  It  is  shown 
theoretically  that  arch  weight  has  significant  effect  on  large  de¬ 
flections.  Results  of  experiments  conducted  are  in  good  agreement  with 


theoretical  results. 
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UN 


LN 


UC 


Young's  modulus 

force  at  the  point  of  loading 

=  F  R  2/EI 
0  0 

rise  of  circular  arch 

moment  of  inertia 

bending  moment  at  any  point 

bending  moment  at  the  point  of  loading 

=  MR  /El 
o  o 

normal  force  at  the  point  of  loading 
applied  concentrated  radial  load 
upper  buckling  load 
lower  buckling  load 

2 

dimensionless  upper  buckling  load  =  P^Ro  /E I 

2 

dimensionless  lower  buckling  load  =  P^Rq  /El 
dimensionless  upper  buckling  load  considering  effects 
of  arch  weight 

dimensionless  lower  buckling  load  considering  effects 
of  arch  weight 

shearing  force  at  the  point  of  loading 
radius  of  curvature  in  the  distorted  state 
radius  of  circular  arch 
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=  s/Rq 

thickness  of  arch 
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PL 
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angle  between  X  axis  and  the  line  of  action  of  F 
angle  between  Y  axis  and  the  normal  at  S 
4>  in  the  undistorted  state 
semi-angle  of  circular  arch 
central  deflection 

-  «/Ro 

weight  per  unit  length  of  the  arch 

-  pr03/ei 

=  p/^UN 

=  ^Uc/^UN 

=  ^LC^UN 

=  R  ft2/t 
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Other  symbols  are  defined  when  used. 
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CHAPTER  I 


INTRODUCTION 

The  shallow  arch  represents  one  of  the  simplest  realistic 
structures  from  which  most  of  the  features  of  elastic  instability  theory 
can  be  illustrated.  The  study  of  this  subject  gives  a  better  insight 
to  the  more  complicated  problem  of  buckling  of  shells  which  has  pre¬ 
sented  a  challenge  to  engineers  for  many  years. 

1.1  Definition  of  Buckling  Load 

Instability  in  the  case  of  shallow  arches  is  caused  by  the 
load  reaching  a  maximum  as  illustrated  in  Figure  1.1,  a  typical  load- 
deflection  curve  for  a  shallow  arch.  The  load-deflection  curve  implies 
that  there  exist  three  possible  states  of  equilibrium  for  each  load  in 
the  range  \  <  ^  <  PU  and  on^  one  eclu^ ibrium  state  for  each  load  out¬ 
side  this  interval.  The  points  on  the  branch  OA  of  the  curve  correspond 
to  unbuckled  states,  those  on  the  branch  BC  to  buckled  states,  and  those 
on  the  segment  AB  to  unstable  states.  The  arch  buckles  at  some  load  in 
the  interval  Py  £  P  £  Py  and  hence  Py  and  Py  are  called  the  lower  and 
upper  buckling  loads  respectively. 

1.2  Review  of  Related  Work 

Buckling  of  shallow  arches  under  concentrated  radial  load  has 
been  investigated  by  two  different  non-linear  methods. 


1 


2 


FIGURE  1.1  A  TYPICAL  LOAD-DEFLECTION  CURVE  FOR  A  SHALLOW  ARCH 


1.2.1  Energy  Criterion 

In  1962,  A.  Gjelsvik  and  S.R.  Bodner  (4)*  presented  a  theo¬ 
retical  and  experimental  analysis  of  this  problem.  Their  theoretical 
investigation  was  based  on  the  energy  methods.  The  total  potential 
energy  due  to  a  central  load  was  expressed  in  terms  of  the  tangential 
and  radial  displacement  functions.  The  tangential  displacement  function 
was  eliminated  from  the  energy  expression  employing  the  condition  of 
equilibrium  in  the  tangential  direction,  i.e.,  the  variation  of  total 
potential  energy  with  respect  to  tangential  displacement  must  be  zero. 

An  approximate  series  solution  to  get  the  complete  load-deflection 
curve  and  hence  the  upper  and  lower  buckling  loads  was  also  given.  Thei 
experimental  work  compared  favourably  with  their  theoretical  solution. 
However,  only  arches  having  an  arch  weight  to  upper  buckling  load  ratio 
less  than  0.05  were  tested  experimentally. 

Whereas  Gjelsvik  and  Bodner  considered  the  variation  of  total 
potential  energy  in  tangential  direction  and  then  followed  the  well 
known  Ritz  Method  (10),  H.L.  Schreyer  and  E.F.  Masur  (9)  considered 
the  variation  of  total  potential  energy  both  in  tangential  and  radial 
directions  and  obtained  a  fourth  order  differential  equation  relating 
load  and  deflection.  An  exact  solution  to  this  equation  was  also  given. 
They  considered  both  symmetrical  and  anti symmetri cal  modes  of  defor¬ 
mation  and  found  that  under  concentrated  load  at  the  center,  the  arch 
remains  symmetrical  until  after  upper  buckling  load  P ^  and  though  un- 


*Numbers  in  parentheses  refer  to  references  in  Bibliography. 
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symmetrical  in  the  unstable  region,  returns  to  a  symmetrical  configu¬ 
ration  when  stability  is  again  attained,  as  illustrated  in  Figures 

1.2  -  1.4.  Hence  for  this  problem  the  symmetrical  buckling  criterion 
governs  regardless  of  arch  geometry. 

1.2.2  Elastica  Approach 

J.S.  Kennedy  (6)  derived  a  non-linear  differential  equation 
to  satisfy  the  geometrical  conditions  and  the  equilibrium  state  of  the 
arch  using  the  general  theory  of  the  elastica  (7).  This  was  subsequently 
solved  for  symmetrical  deformations  by  algebraic  transformations  and 
the  introduction  of  elliptic  integrals.  Kennedy  also  conducted  a  series 
of  experiments  with  light,  thin,  circular  arches  which  were  induced  to 
buckle  symmetrical ly  by  clamping  the  centre  of  the  arch.  The  arch 
centre,  therefore,  was  allowed  to  move  only  in  the  vertical  direction 
and  a  horizontal  tangent  at  that  point  was  maintained.  For  a  given 
subtended  angle  of  the  arch  the  application  of  Kennedy's  theory  re¬ 
sulted  in  only  one  load-deflection  curve  for  non-dimensional  plotting 
whereas  he  obtained  several  different,  although  similar,  curves  experi¬ 
mentally  by  employing  various  arch  geometries.  Poor  arch  clamping  was 
suggested  as  the  reason  of  this  discrepancy. 

Later  in  1963,  D.A.  Burns  (2)  undertook  an  experimental  investi¬ 
gation  of  shallow  arches  under  a  concentrated  load.  He  tested  a  number 
of  light  clamped  arches  and  found  that  in  some  cases  there  was  as  much 
as  35%  deviation  from  the  theoretical  curves  obtained  by  the  appli¬ 
cation  of  Kennedy's  or  Bodner's  theory.  The  arch  specimens  tested  by 
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FIGURE  1.2  ARCH  IN  UNBUCKLED  STATE,  RETAINING 

SYMMETRY. 


FIGURE  1.3  ARCH  IN  UNSTABLE  STATE,  SHOWING 

PRONOUNCED  ASYMMETRY. 


FIGURE  1.4  ARCH  IN  BUCKLED  STATE,  SHOWING 

RETURN  TO  SYMMETRY 


6 

Burns  had  weight  to  upper  buckling  load  ratios  from  0.06  to  0.48. 

Arch  weight  was  suggested  as  the  cause  of  the  large  discrepancies  be¬ 
tween  the  theoretical  and  experimental  curves. 

1 .3  Aim  of  the  Thesis 

In  the  present  work,  Burn's  line  of  thought,  i.e,,  weight  of 
the  arch  does  have  an  effect  on  the  buckling  load,  has  been  investigated 
further.  A  non-linear  differential  equation  using  the  elastica  ap¬ 
proach  (7)  and  taking  the  arch  weight  into  consideration  has  been 
derived.  The  equation  developed  has  been  solved  to  obtain  complete 
load-deformation  curves  by  numerical  techniques  described  in  Ap¬ 
pendices  A  and  B.  Experiments  on  arches  covering  a  wide  range  of  arch 
geometry  have  been  conducted  and  the  results  are  in  close  agreement 
with  the  theoretical  curves.  A  linear  correlation  has  been  observed 
between  the  arch  weight  to  upper  buckling  load  ratio  and  non-dimension- 
lised  upper  and  lower  buckling  loads. 
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CHAPTER  II 


THEQRETICAL  ANALYSIS 

The  elastica  approach  employed  to  study  the  large  deflections 
of  flexible  bars  is  well  known  (7).  This  method  is  used  here  to 
analyse  the  symmetrical  mode  of  deformation  of  a  ring  segment  of  radius 
RQ,  having  built  in  ends  which  subtend  an  angle  2ft,  loaded  by  a  radial 
load  P  directed  along  the  line  of  symmetry  as  shown  in  Figure  2.1.  The 
plane  of  loading,  i.e.,  the  plane  defined  by  the  applied  load  and  the 
reactive  forces,  contains  one  principal  axis  of  inertia  of  each  section, 
so  that  the  bending  of  the  arch  is  confined  to  that  plane  only. 

2.1  Derivation  of  Governing  Equation 

Because  of  the  symmetry,  only  one  half  of  the  arch  need  be 
considered  in  the  development  of  the  governing  differential  equation 
for  the  problem.  As  shown  in  Figure  2.2,  at  the  mid  point  SQ,  the  load 
will  consist  of  a  couple  M  ,  a  horizontal  compressive  force  NQ  and  a 
vertical  shearing  force  Q  ,  and  at  some  other  point  S  of  such  a  couple 
and  a  force  as  to  provide  equilibrium.  A  coordinate  system  X,Y  with 
its  origin  at  the  load  point  S  ,  the  X  axis  coinciding  with  the  tangent 
at  Sq  and  the  Y  axis  coinciding  with  the  normal  at  SQ  and  lying  along 
the  radius  of  curvature  will  be  employed.  The  forces  Nq  and  Q0  acting 
at  the  point  SQ  can  be  replaced  by  a  single  force  FQ  acting  along  a 
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FIGURE  2.1  CLAMPED  CIRCULAR  ARCH  UNDER  CONCENTRATED  LOAD 
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FIGURE  2.2  SEMI-ARCH  AFTER  DEFORMATION 
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line  which  makes  an  angle  3  with  the  X  axis. 

The  bending  moment  at  an  arbitrary  point  is  defined  as  positive 
if  it  tends  to  decrease  the  curvature.  From  equilibrium  the  moment  M, 
at  any  point  S,  on  the  curve  is  given  by 


M  =  Mq  -  FQ(xsin3  -  ycos3) 


x_.=x 


x.j=0 


p(x-x.) 
cos<|> . 


2.1 


where  p  is  the  weight  per  unit  length  of  the  arch  and  <j>.  is  the  angle 
between  the  Y  axis  and  the  normal  at  . 

Also,  if  R  is  the  radius  of  curvature  in  the  distorted  state 
at  the  arbitrary  point  S,  the  moment  at  that  point  is 


M  =  EI(J-  -  1)  2.2 

where  El  is  flexural  rigidity  of  the  arch. 

Passing  from  the  X,Y  coordinates  to  the  natural  coordinates 
s  and  <j>,  where  s  is  the  length  of  the  arch  SqS,  and  <f>  is  the  angle  be¬ 
tween  the  Y  axis  and  the  normal  at  S,  and  letting  <|>  =  <J>  in  the  un¬ 
distorted  state,  there  exist  the  relations 


1  =  d£ 

R  ds 


d x 
ds 


COScj) 


sin* 


2.3 


ar  jrf’oa  tsrtt  3  b  3ns  10m  ,2  JnFoq  9ri$  Jfi 


.  n  9ri  1  t-o  r  i  r  '&  U  Mt  ai  13  y.srti 


■  ulgf  j  13  6  ,22  fiD*>  F  o  rttgnsf  .  .r  j  v  ,  w  «4>  bns  a 

4>  pnrtJsf  bfi  >  ,2  te  f  sun  on  9iiJ  b ns  arxf  Y  9riJ  n?9wJ 
anofj  ft  i  sriJ  Jcrxe  atari/  t9^6Ja  D9/to:Lrb 


11 


Employing  Equations  2.2  and  2.3;  Equation  2.1  becomes 


d<f> 


EI  (35s  -  ^)  =  -  Fjxsing  -  ycosB) 


ds 


X.j  =  X 


x^O 


p(x-X.) 

cos4~ 


dx.j 


2.4 


The  slope  of  the  tangent  cj) . ,  is  a  complicated  function  of  the 
loading.  Therefore,  the  last  term  in  Equation  2.4  which  accounts  for 
the  direct  contribution  of  the  arch  weight  to  the  bending  moment  becomes 
difficult  if  not  impossible  to  evaluate.  To  simplify  the  situation,  let 
us  restrict  our  attention  to  shallow  arches,  i.e.,  arches  in  which  the 
ratio  of  height  H,  to  radius  of  curvature  R  ,  is  less  than  0.1.  For 
such  cases,  in  evaluating  the  last  term  of  Equation  2.4,  it  is  reasonable 
to  employ  the  approximation  that  each  point  on  the  arch  moves  vertically. 
This  means  that  the  value  of  the  last  term  of  Equation  2.4  at  a  section 
will  always  be  the  same  regardless  of  the  arch  deformation.  Consequently, 
the  bending  moment  due  to  the  arch  weight  can  be  evaluated  with  respect 
to  the  undistorted  state. 

<v 

Let  X,  Y  be  the  coordinate  axes  fixed  to  the  undistorted  shape 

of  the  arch  and  d>  be  the  angle  between  the  Y  axis  and  the  normal  to  S, 

0 

as  shown  in  Figure  2.3.  Then  the  last  term  in  the  Equation  2.4  can  be 


written  as 


.gntbfior 
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FIGURE  2.3  SEMI-ARCH  BEFORE  DEFORMATION 


5=<t> 

o 

/  pR  (x  -  R  sin£)  d£ 
5=0  0  0 


5=cj> 

=  [pRn(x5  +  R  cosC)]  0 
u  u  5=o 


=  pr0(M>0  -  y) 


and  the  Equation  2.4  now  takes  the  form 


dch 


EI^ds^  "  Hs)  =  Mo  “  F0(xsin13  “  ycos3) 


-  pr0(M>0  -  y) 


2.5 


Differentiating  with  respect  to  s,  the  Equation  2.5  becomes 


EI(^'$}  =  '  Fo(^SinS  -^C0se) 


~  d4 

- ;  ar  -  f ' 


Employing  the  relations  2.3  and  noting  that  the  curvature  of  ring  seg¬ 
ment  is  uniform  in  the  unloaded  state,  i.e.,  d(f)o/ds  =  1  / Rq ,  Equation 
2.6  simplifies  to 


=  FQ  sin(<j>-$)  -  pRQ(f)ocos4)0 
ds 


-  El 


2 . 7 


' 


)u«  .  n  lo  9' u:  'uid  ■;*  Jfcrtf  gn  r  torr  brtf>  anotJr>r#ni  9fJ  pn  f  qni3 


o  upJ  r  fl\  :j\ib  <.s  r  1 9  oJ  z  t  r>r,  c  iu  ®rto  rr  mroi  hu  2f  Jnsm 
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Adopting  Rq  and  El  as  standards  of  length  and  stiffness  and 
defining  the  dimensionless  quantities, 


MR 
o  o 

El 


F  = 


FoRo‘ 

El 


3 


El 


Equation  2.7  can  be  written  in  the  form 


+  Fs i n ( )  =  p  s  cos  s  2.8 

ds2 

Equation  2.8  is  the  differental  equation  describing  load- 
deformation  characteristics  of  the  system.  It  may  be  remarked  here 
that  this  equation  is  of  the  same  form  as  developed  by  Kennedy  (6)  but 
in  addition  it  has  an  extra  term  on  the  right  hand  side  which  accounts 
for  the  effects  of  arch  weight. 

2.2  Boundary  Conditions 

Since  we  are  analysing  the  symmetrical  mode  of  deformation 
of  the  arch,  the  tangent  at  the  center  must  remain  horizontal,  i.e., 


cj)  =  0  at  s  =  0 


2.9 


sett  aarJan  Jos^rb  flottsmcrtsb 
iu  1  (d  J  Yb9fir.j>  Yd  b'-qo  svafc  25  mio^t  sn^a  *rtt  ^o  t\  riof*  ups  ai  IJ  3btii 
a  “  '-5  bl  3  -  I  ns  'no  fn'  -  r  i  K«  i  c  5  j'r  fi  i  n  nf 
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Secondly,  because  of  the  fixed  ends 

<J>  =  ft  at  7  =  ft  2.10 

In  addition  to  the  above  two  boundary  conditions,  the  hori¬ 
zontal  projection  must  remain  constant  and  equal  to  sin  ft,  i.e., 


ft  _ 

/  —  ds  =  sin  ft 

0  ds 

Since 

~  =  cos  <j>  , 
ds 

then 

ft 

/  cos^  ds  =  sin  ft  2.11 

0 

A  solution  to  the  governing  differential  equation  2.8  sub¬ 
ject  to  the  boundary  conditions  2.9  -  2.11  is  desired.  Once  a  solution 
has  been  obtained,  the  non-dimensional  central  deflection  6,  at  the 
point  of  loading  is  given  by 


Since 

then 

ft  i — 

Ro  0  ds 

— ^  =  sin  (p  , 
ds 

ft 

6  =  (1  -  cos  ft)  -  /  sin<j>  ds 

0 

then 


:  CO  -fc  -  )d  3\  :  oi  vo  : 


930  l 


-c!j2  8.S  norJsups  Cfif:J,n9't9^ fb  gnfrfmvog  9fid  ctf  nofJufoa  A 

. 


.  f  •  •  : 
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The  non-dimensional  concentrated  load  parameter  is  given  by 


P  =  2F  sin  3 


2.13 


2.3  Deflection  Under  Dead  Weight 

The  central  deflection  because  of  the  arch  dead  weight  was 
found  by  application  of  the  linear  theory  in  form  of  Castigl iano's 
Theorem. 

Considering  one  half  of  the  arch,  Figure  2.3,  the  moment  at 
any  point  S,  is  given  by 


M  =  M  +  N  R  ( 1  -cos4)  )  -  Q  R  sin  <b 
0  00  T  o  ^00  ”0 


2.14 


5=0 


which  on  simplification  gives 


2.15 


By  Castigl iano 's  Theorem 


2.16 


ro^oarir  2'on6rr9rj260 
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R  £1  r.M 

%  =  ITl0"Wo**o 

are  respectively  the  rotation,  the  horizontal  deflection  and  the  vertical 
deflection  of  the  center  point. 

Substituting  Equation  2.15  into  the  Equations  2.16,  carrying 
out  the  required  integration  and  taking  the  limit  of  the  resulting  ex¬ 
pressions  as  Qq  approaches  zero  yields, 

R 

6Mo  =  ET  +  NoRo(fi  '  sinn) 

-  pR  2(2sinfi!  -  ficosfi  -  fl)]  2.17 

R  2 

6N0  =  [Mo(fi-sinfi)  +  NoRq(|^  +  -  2sinft) 

p  3  3  1 

-  pR  (3sin£2  -  Slcosfl  -  j  SJ  -  g  sin2f2  +  flcos2fl)]  2.18 

R  2 

6Q0  =  [M0(cosn-D  +  NoRo(cosn  -  |  -  ^2|2£) 

+  pR  2(^-  -  ^  S7sin2fi  -  |  cos2fi  +  cosfl  -  |)]  2.19 

The  symmetry  of  the  problem  requires  that 

6RI  =  6N  =  0 
o  o 


2.20 


(finfaS  -  -  +  ;  +  (W-9)aM]  4;  * 


18 


Therefore  from  Equations  2.17  and  2.18,  Mq  and  NQ  can  be  found 
in  terms  of  p,  ft,  Rq  and  El.  Knowing  the  values  of  MQ  and  N  ,  6Qq,  the 
central  deflection  under  dead  weight,  can  be  calculated.  For  the  arches 
considered  in  this  thesis,  this  deflection  was  found  to  be  negligible, 
i.e.,  less  than  one  thousandth  of  an  inch. 


CHAPTER  III 


METHOD  OF  SOLUTION 


3.1  Integration  Procedure 

The  governing  differential  equation,  2.8,  can  be  written  as 
a  system  of  two  simultaneous,  first-order  differential  equations 

dd) 

-f;  =  a 
ds 

~  =  ps  cos  s’  -  T  sin(cj>-$)  3.1 

ds 

Such  systems  can  usually  be  solved  numerically  without  diffi¬ 
culty,  inspite  of  non-linear  terms,  provided  that  the  problem  is  of  the 
11  initial -value  type",  with  cf>  and  a  known  at  some  "initial"  value  of  ?. 

In  the  present  case,  at  the  initial  point  F  =  0,  <j>  is  zero  but  a  and  the 
angle  3,  which  is  dependent  on  the  loading  are  unknown.  To  overcome 
this  problem,  the  author  has  used  a  "Shooting  Method".  The  slope  a  and 
the  angle  3  are  assumed  at  the  initial  point,  Equations  3.1  are  integrated 
numerically  using  Gill's  modification  of  the  Runge-Kutta  fourth  order 
method  (Appendix  A)  and  the  terminal  conditions  2.10  and  2.11  computed. 

The  assumed  initial  values  of  a  and  3  are  then  adjusted  iteratively  until 
the  boundary  conditions  are  satisfied. 
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3.2  Block  Diagram 

The  computation  procedure  is  illustrated  in  the  block  diagram. 
Figure  3.1.  For  convenience,  the  following  notation  will  be  employed. 

<f>  =  4>  at  ?  =  0 
o 

4>  =  4>o  s’  =  ft 

a  =  a  at  s  =  0  3.2 

o 

The  computation  scheme  incorporates  two  iterative  loops  be¬ 
cause  of  the  two  unknown  quantities  3  and  a  .  Corresponding  to  an 
assumed  value  of  angle  3,  let  a-j ,  be  two  guesses  of  the  initial  a  , 
and  let  ^(a-j),  ^(c^)  be  two  values  of  <j>  at  7  =  ft  obtained  from  inte¬ 
grating  the  differential  equations  3.1.  Graphically  the  situation  may 
be  presented  as  in  Figures  3.2  and  3.3.  In  Figure  3.2,  the  solutions 
of  the  initial-value  problem  are  drawn,  while  in  Figure  3.3,  <|>^(a  )  is 
plotted  as  a  function  of  a  .  Normally  a  better  approximation  to  aQ 
can  now  be  obtained  by  linear  interpolation.  The  intersection  of  the 
line  joining  P-j  and  with  a  horizontal  line  from  <|>^(a  )  =  ft,  the  re¬ 
quired  value  at  s  =  ft,  has  as  its  aQ  coordinate 

ou  =  an  +  (a0-an )  t— ? — ? - t — ? — r 

3  1  2  i  vFw 


3.3 


r  Qx  ti  l  ^u|  M  m  9f  riw  .riwb-ib  e'ts  mafdoiq  suFfiv-rbrd  i  nr  ertt 
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FIGURE  3.1  BLOCK  DIAGRAM 


FIGURE  3.2  <p  vs  s 


FIGURE  3.3  4^(01^)  vs  % 
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The  system  of  differential  equations  3.1  is  now  again  integrated 
using  the  initial  value  aQ  =  to  obtain  ^(c^).  Again  using  linear 
interpolation  based  on  a^,  a^,  a  new  approximate  is  obtained.  The 
process  is  repeated  until  the  convergence  has  been  obtained  i.e.  the 
boundary  condition  2.10  is  satisfied  to  the  desired  accuracy.  In  a 
similar  manner,  the  angle  3  is  corrected  to  satisfy  the  condition  2.11. 

The  rapidity  of  convergence  is  clearly  dependent  upon  how  close  the 
initial  guesses  are  to  the  actual  values.  It  may  be  noted  here  that  any 
change  in  the  angle  3,  will  change  the  slope  a  ,  therefore  the  inner 
iterative  loop  will  have  to  be  satisfied  again. 

The  central  deflection  6,  given  by  the  equation  2.12  is  then 
computed  by  stepwise  integration  employing  a  combination  of  Simpson's 
and  Newton's  Three-Eighth  Rule  (Appendix  B). 
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CHAPTER  IV 


EXPERIMENTAL  PROCEDURE  AND  APPARATUS 

4.1  Apparatus 

The  experimental  apparatus  is  shown  in  Figure  4.1. 

The  arch  bed  was  machined  from  a  length  of  65S-T6  6"  x  2-3/4" 
aluminum  channel.  The  top  face  and  the  tips  of  channel  flanges  were 
milled  in  the  Department  Workshop  to  ensure  that  these  surfaces  would 
be  smooth  and  parallel  to  each  other.  A  centre  slot  was  machined  out 
of  the  top  face  of  the  bed  as  shown  in  Figure  4.2  to  permit  large  de¬ 
formations  of  the  arch  should  this  prove  necessary  during  a  test.  The 
parallel  rows  of  holes  were  drilled  so  that  the  arch  support  blocks 
would  be  positioned  via  guide  pegs,  and  would  thus  be  aligned  correctly 
with  no  further  effort.  The  holes  tapped  at  regular  intervals  along 
the  channel  centerline  were  placed  to  mate  with  the  fixing  bolt  of  each 
support  block. 

Three  sets  of  support  blocks  were  fabricated.  These  blocks 
were  milled  so  that  when  the  clamps  were  set  in  position  on  the  arch 
bed,  the  slope  of  the  built-in  end  of  the  arch  would  be  10°,  15°  and 
20°  to  the  horizontal.  Details  of  the  support  block  are  depicted  in 
Figure  4.3. 

The  arch  bed  was  fixed  to  the  scale  pan  of  an  Ohaus  Triple  Beam 
Balance.  The  balance  had  been  securely  located  on  a  firm  wooden  bench 
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FIGURE  4.1  GENERAL  VIEW  OF  APPARATUS. 


1/4"  GUIDE  HOLES  CENTERAL  SLOT  1/4"  TAPPED  HOLES 
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FIGURE  4.3  SUPPORT  BLOCK 
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and  shimmed  until  level.  The  bed  scale  system  was  adjusted  carefully 
and  it  was  checked  with  a  sensitive  spirit  level  to  ensure  that  the 
upper  surface  of  the  arch  bed  lay  in  the  horizontal  plane. 

Arch  deflections  were  applied  by  a  loading  bar  fashioned 
from  a  hook  gauge,  which  was  supported  by  a  tubular  steel  frame.  A 
Vernier,  with  least  count  0.012  inch  was  used  with  the  loading  bar. 

The  central  point  of  the  arch  was  allowed  to  move  only  in  the  vertical 
direction  and  the  tangent  at  the  midpoint  was  maintained  horizontal 
by  the  arrangement  shown  in  Figures  4.4  and  4.5. 

Prior  to  its  installation,  the  Ohaus  Triple  Beam  Balance  was 
calibrated,  using  weights  supplied  with  an  analytical  balance,  accurate 
to  a  milligram.  It  was  found  to  be  within  the  specified  tolerance  of 
one  gram. 

4.2  Test  Specimens 

The  arch  test-pieces  were  cut  from  banding  steel  0.015  or 
0.026  inch  thick  and  0.75  inch  wide.  The  required  length  was  calculated 
from  the  specified  angle  and  radius  of  curvature  of  the  arch.  Centre¬ 
lines  were  scribed  on  each  piece  to  correspond  to  similar  marks  on  the 
supporting  blocks  thereby  assisting  in  the  alignment.  The  black  finish 
of  the  banding  steel  lent  itself  readily  to  this  process.  Each  length 
was  carefully  inspected  for  imperfections  and  kinks  prior  to  testing. 

The  circularity  of  the  arch  was  further  insured  by  measuring  the  rise 
of  the  arch  by  a  Vernier  Height  Gauge  accurate  to  one  thousandth  of  an 


inch. 
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FIGURE  4.4  CENTRAL  CLAMP. 


FIGURE  4.5  DEFORMED  ARCH  RETAINING  SYMMETRY. 
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Pertinent  data  for  the  eleven  arches  which  were  tested  experi¬ 
mentally  is  given  in  Table  4.1,  and,  for  the  sake  of  convenience,  each 
specimen  has  been  given  an  identifying  number. 

The  arches  selected  covered  a  wide  range  of  arch  geometry  and 
weight  to  upper  buckling  load  ratio. 

4.3  Experimental  Procedure 

With  a  test-piece  firmly  located  in  the  apparatus,  the  weight 
required  to  balance  the  scale  was  noted.  The  loading  bar  was  lowered 
until  the  loading  head  just  touched  the  arch.  The  head  was  securely 
fixed  to  the  arch  center  to  keep  the  tangent  at  that  point  horizontal. 
The  arch  was  then  deflected  incrementally,  and  for  each  setting  the 
balancing  weight  registered  on  the  scale  arm  was  noted.  Thus  with  the 
scale  adjusted  so  that  the  scale  pan  was  again  held  at  the  original 
height,  a  situation  was  obtained  whereby  a  measured  central  load  in¬ 
duced  the  arch  to  deflect  a  given  distance.  The  zero  reading,  i.e.,  the 
weight  required  to  balance  without  any  deflection  was  again  noted,  and 
if  it  differed  from  the  initial  value,  the  whole  set  of  readings  was 
discarded.  The  process  was  repeated  three  to  four  times  for  each  speci¬ 
men  and  the  average  readings  were  plotted. 


;  ;■>  ii  ■  ' .  S’  :j  ' 


>g  £6  . .  ftsb  rU  'rw  sans  fed  o3  bsifupsT  drip  tew 


c  o  >  9"  9^6'i  tv  \  add  .j ns  n9m 


31 


Table  4.1 

Dimensions  of  Test  Specimens 


Identifying 

Number 

Semi -angle 

Radius  of 

Curvature 

R„  in  inches 

0 

Thi ckness 
t  in 
inches 

Non-dimensional 

parameter 

X  =  RQf22/t 

SI 

10° 

82.67 

0.015 

166.77 

S2 

10° 

69.08 

0.015 

139.36 

S3 

10° 

82.67 

0.026 

96.86 

S4 

10° 

69.08 

0.026 

80.94 

S5 

15° 

55.47 

0.015 

251.76 

S6 

15° 

46.35 

0.015 

210.38 

S7 

15° 

55.47 

0.026 

146.22 

S8 

15° 

46.35 

0.026 

122.18 

S9 

o 

O 

C\J 

41.97 

0.015 

338.70 

S10 

20° 

41  .97 

0.026 

196.71 

Sll 

o 

O 

C\J 

35.07 

0.026 

164.37 

Width  =  0.75  inch 

Modulus  of  Elasticity  =  29.9  x  10^  psi 
Material:  Banding  Steel 
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CHAPTER  V 


RESULTS  AND  CONCLUSIONS 

5.1  Load-Deflection  Curves 

The  load-deflection  curves  for  the  symmetrical  deformation  of 
specimens  SI  -  Sll,  obtained  by  the  application  of  the  elastica  approach 
neglecting  arch  weight,  i.e.,  by  the  solution  of  differential  equation 
as  outlined  in  (1),  as  well  as  the  ones  obtained  by  considering  the 
effects  of  arch  weight,  i.e.,  by  the  solution  of  differential  equation 
2.8,  are  presented  in  Figures  5.1  to  5.11.  Experimental  results  ob¬ 
tained  by  the  deflection  controlled  procedure  explained  in  Chapter  IV 
are  included  for  each  specimen.  In  all  cases,  results  obtained  from 
the  Gjelsvik  and  Bodner  (3)  and  the  Schreyer  and  Masur  (8)  theories  are 
given  for  comparison.  Whereas  the  experimental  results  have  a  large 
discrepancy  in  most  of  the  cases  from  the  previous  work  of  Gjelsvik  and 
Bodner,  Schreyer  and  Masur,  and  Kennedy,  they  are  in  good  agreement 
with  the  present  theoretical  analysis. 

5.2  Buckling  Loads 

Upper  and  Lower  buckling  loads  for  specimens  SI  -  Sll  obtained 
from  the  respective  Load-Deflection  curve  are  given  in  Table  5.1. 
Following  notations  are  used 

=  Non-dimensional  Upper  Buckling  Load 


32 


' 


i  D6oviq?>  eor^&fs  ertj  rtor^eohq  6  (c  b*>'  6Jdo  .  ff 2  -  f2  ansmfosqa 


!  •  !  > 


' 


;.)Q2  '  t  ,  of  gnobuo  •v?v."J  „n6  9Q  U 


t  u  Ji6  anc  J63on  gniwoffol 


33 


FIGURE  5.1  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  SI 
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FIGURE  5.2  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S2 
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FIGURE  5.3  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S3 
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FIGURE  5.4  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S4 
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FIGURE  5.5  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S5 
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Neglecting  Weight 
Experimental 
Gjelsvik  &  Bodner's 
Schreyer  &  Masur's 


Theory 
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FIGURE  5.6  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S6 
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FIGURE  5.7  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S7 
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FIGURE  5.8  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S8 


seo.o 
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FIGURE  5.9  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S9 
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FIGURE  5.10  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  S10 
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FIGURE  5.11  LOAD-DEFLECTION  CURVES  FOR  SPECIMEN  Sll 
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SUMMARY  OF  THEORETICAL  RESULTS 
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PLN  =  Non_dl'mens''ona'l  Lower  Buckling  Load 
PUC  =  Non-dimensional  Upper  Buckling  Load  considering  the 
effects  of  arch  weight 

P|_q  =  Non-dimensional  Lower  Buckling  Load  considering  the 
effects  of  arch  weight 

Pw  =  Ratio  of  non-dimensional  arch  weight  and  Upper  Buckling 
Load 


pu 

pL 


p/PUN 

=  Ratio  of  P”uc  and 


PUC/PUN 


=  Ratio  of  and  P ^ 


PLC/PUN 

py  and  are  plotted  as  functions  of  p^  in  Figure  5,12. 


5,3  Conclusions 

1,  Figures  5,1  -  5.11  show  that  the  arch  weight  has  a  significant 
effect  on  the  large  deflections  of  thin  arches  under  concentrated  load. 

2,  There  is  a  linear  relationship  between  p  and  pn  and  between 

w  u 

Pw  and  p^,  as  shown  in  Figure  5,12,  Both  upper  and  lower  buckling  loads 
decrease  as  the  arch  weight  to  upper  buckling  load  ratio  increases. 

3,  Keeping  the  semi-angle  and  stiffness  of  the  arch  constant,  the 
arch  weight  to  buckling  load  ratio  will  increase  by  increasing  the  span. 
Therefore,  in  the  case  of  longer  arches,  arch  weight  will  have  a  more 
significant  effect  on  the  buckling  loads. 


' 
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FIGURE  5.12  p|j  AND  p.  vs  pw 
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4.  For  specimens  SI,  S2,  S5  and  S9,  the  lower  buckling  load  is 
lowered  to  the  extent  that  it  becomes  less  than  zero.  Therefore  it  is 
possible  for  these  arches  to  remain  in  a  buckled  configuration  under 
their  own  weight.  This  phenomenon  was  also  observed  experimental ly . 

On  being  given  a  small  disturbance,  specimens  SI,  S2,  S5  and  S9  jumped 
to  the  symmetrical  buckled  shape  and  stayed  in  that  form,  Figure  5.13, 
even  after  the  load  causing  the  disturbance  had  been  removed. 

For  any  arch  having  arch  weight  to  upper  buckling  load  ratio 
of  more  than  0.53,  Figure  5,12,  the  lower  buckling  load  is  less  than 
zero.  Therefore  this  phenomenon  of  arch  buckling  under  its  own  weight 
when  subjected  to  a  disturbance  is  possible  for  any  arch  having  a 
weight  to  upper  buckling  load  ratio  of  more  than  0.53, 

5.  Gjelsvik  and  Bodner  (4)  were  successful  in  obtaining  a  good 
agreement  between  their  experimental  and  theoretical  results  because 
for  the  arches  tested  by  them  the  arch  weight  to  upper  buckling  load 
ratio  was  small,  i.e.,  less  than  0.05  and  hence  the  assumption  of 
negligible  arch  weight  was  permissible. 
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FIGURE  5.13  BUCKLED  ARCH  UNDER  ITS  OWN  WEIGHT. 


CHAPTER  VI 


FURTHER  APPLICATIONS  OF  THE  "SHOOTING"  METHOD 


1.  The  "shooting"  method  used  by  the  author  presents  an  elegant 
tool  for  the  study  of  the  large  deflections  of  flexible  elastic  bars. 

A  number  of  non-linear  boundary  value  problems  can  be  solved  by  this 
technique.  For  example,  consider  the  large  deflections  of  complete 
circular  ring  subjected  to  two  equal  and  opposite,  concentrated  loads 
acting  along  a  diameter.  The  complete  ring  is  free  to  expand  in  the 
direction  at  right  angles  to  the  line  of  action  of  the  loads.  There¬ 
fore  there  will  be  no  thrust  N  ,  Figure  2.2.  The  angle  3  will  there¬ 
fore  be  tt/2.  The  Equation  2.8  for  this  case  reduces  to 


ds 


6.1 


neglecting  effects  due  to  ring  weight. 

Boundary  conditions  for  this  case  are 


0=0  at  s  =  0 


6.2 
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In  this  case  the  unknown  initial  value  is  only.  Follow- 

ds 

ing  the  same  procedure  as  before  the  complete  load  deflection  curve 
is  obtained,  Figure  6.1,  For  comparison,  results  obtained  by  Timeshenko 
and  Gere  (10)  are  also  given  in  Figure  6.1, 

2.  With  suitable  modifications  in  boundary  conditions,  the 
method  used  in  this  thesis  can  also  be  applied  to  solve  symmetrical 
deflections  of  pinned  ends  arches  under  concentrated  radial  load  at 
the  point  of  symmetry.  Conditions  2.9  and  2.11  remain  unchanged. 

Since  there  is  no  moment  at  the  pinned  ends,  boundary  condition  2,10 
must  be  replaced  by 


=  ^2.  =  1  at  s  =  n  6.3 
ds  ds 

3.  Arches  considered  in  this  thesis  were  only  0.75  inch  wide. 
The  present  work  can  be  extended  to  cover  wider  arches  and  possibly 
cylindrical  shells  taking  anti  clastic  deformation  into  account. 
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APPENDIX  A 


NUMERICAL  SOLUTION  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 
A.l  Introduction 

The  system  of  differential  equations,  3.1,  has  been  solved 
by  Gill's  variation  of  the  Runge-Kutta  Fourth  Order  method  (8).  This 
is  a  self  starting  integration  procedure  which  is  applicable  to  a 
system  of  n-first  order  ordinary  differential  equations  and  is  highly 
suitable  for  use  on  a  digital  computer. 

A. 2  Runge-Kutta  Methods 

Given  the  system  of  first-order  ordinary  differential  equations 

dy  • 

=  y.j  =  f^x.y-jU),  y2(x),  ...  yn(x) ) 
with  the  initial  conditions 


y •  (x  )  =  y . 

J  i  o  J  i  o 

it  is  desired  to  obtain  the  values 

y  i  ( x0  +  h) 
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where  h  is  an  increment  of  the  independent  variable  x. 

The  basis  of  all  Runge-Kutta  Methods  is  to  express  the  differ¬ 
ence  between  the  values  of  y  at  xn+1  and  xn  as 

m 

Vi  -  yn  =  il1ciki 

where  c.  are  constants  and 

i-1 

ki  =  hn  f(xn  +  aihn>  +  .1  6ijkj>  • 

J  , 

c^ ,  a..  and  3^  are  obtained  to  have  desired  properties  by 
approximating  the  Taylor  series  solutions. 

A.  3  Gill's  Modification 

To  apply,  the  method  on  high-speed  digital  computer.  Gill  (3) 
developed  a  calculation  procedure  in  which  the  constants  c^ ,  a.,  and 

B.  .  are  obtained  to  minimize  the  storage  problem  and  to  give  the  highest 

*  J 

attainable  accuracy.  In  addition.  Gill's  method  requires  comparatively 
few  instructions  and  hence  results  in  substantial  savings  in  computer 
time.  The  calculation  proceeds  by  sequentially  calculating  the  j  ,  k  , 
and  q  values  as  given  in  Table  A.l. 

A. 4  Accuracy  of  the  Method 

All  the  calculations  are  carried  out  in  double-precision, 
that  is  carrying  16  significant  figures.  Therefore  it  is  reasonable 
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Runge-Kutta  Gill  Scheme  for  First  Order  Differential  Equation 
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to  assume  that  round  off  errors  are  negligible. 

The  truncation  errors  in  the  process  are  caused  because  of 
neglecting  terms  0(h5)  and  higher  in  the  Taylor  Series  expansion  about 
the  initial  values.  The  evaluations  of  higher  order  terms  is  tedious. 

So  it  is  not  possible  to  estimate  per  step  error.  In  order  to  over¬ 
come  this  disadvantage,  control  of  accuracy  is  accomplished  by  adjusting 
the  increment  in  x.  A  comparison  is  made  between  the  function  values 
obtained  at  x  =  xq  +  2h  by  firstly  using  the  increment  h  in  two  stages 
and  secondly  by  using  double  the  increment  2h  in  a  single  stage. 

Let: 

5 

c-j  h  denote  the  truncation  error  using  step  h 
c2(2h)  denote  the  truncation  error  using  step  2h 

y^  denote  the  value  obtained  at  xQ  +  2h  when  using  step  h 
y^2^  denote  the  value  obtained  at  xq  +  2h  when  using  step  2h 
y  denote  the  true  value  of  y  at  xQ  +  2h  . 

Then 

y  -  y^  =  2  Ci  h5 
y  -  y^  =  c2( 2h )5 
For  small  h,  c-|  *  c2,  therefore 

y  -  y(1)  =  js  [y(2)  -  y(1)] 

which  gives  a  measure  of  the  error  introduced.  To  keep  a  check  on  the 
error,  this  procedure  has  been  included  in  the  integration  procedure. 
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is  compared  to  the  predefined  accuracy 


The  introduced  error  y  -  y^ 

(10"  in  the  present  problem)  and  if  greater  than  desired,  the  step 
size  h,  is  decreased.  A  Fortran  IV  subroutine  DRKGS,  based  on  the 
method  described  above  is  given  in  section  A. 5. 
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A. 5  SUBROUTINE  DRKGS 
PURPOSE 

TO  SOLVE  A  SYSTEM  OF  FIRST  ORDER  ORDINARY  DIFFERENTIAL 
EQUATIONS  WITH  GIVEN  INITIAL  VALUES. 


USAGE 

CALL  DRKGS  (PRMT ,Y  ,DERY ,NDIM, IHLF,FCT,OUTP ,AUX) 
PARAMETERS  FCT  AND  OUTP  REQUIRE  AN  EXTERNAL  STATEMENT. 


DESCRIPTION  OF  PARAMETERS 


PRMT 


PRMT ( 1 ) 
PRMT (2) 
PRMT(3) 

PRMT (4) 


PRMT (5) 


Y 


DERY 


DOUBLE  PRECISION  INPUT  AND  OUTPUT  VECTOR  WITH 
DIMENSION  GREATER  THAN  OR  EQUAL  TO  5,  WHICH 
SPECIFIES  THE  PARAMETERS  OF  THE  INTERVAL  AND  OF 
ACCURACY  AND  WHICH  SERVES  FOR  COMMUNICATION  BETWEEN 
OUTPUT  SUBROUTINE  (FURNISHED  BY  THE  USER)  AND 
SUBROUTINE  DRKGS.  EXCEPT  PRMT(5)  THE  COMPONENTS 
ARE  NOT  DESTROYED  BY  SUBROUTINE  DRKGS  AND  THEY  ARE 
LOWER  BOUND  OF  THE  INTERVAL  (INPUT), 

UPPER  BOUND  OF  THE  INTERVAL  (INPUT), 

INITIAL  INCREMENT  OF  THE  INDEPENDENT  VARIABLE 
(INPUT), 

UPPER  ERROR  BOUND  (INPUT).  IF  ABSOLUTE  ERROR  IS 
GREATER  THAN  PRMT(4),  INCREMENT  GETS  HALVED. 

IF  INCREMENT  IS  LESS  THAN  PRMT(3)  AND  ABSOLUTE 
ERROR  LESS  THAN  PRMT(4)/50,  INCREMENT  GETS  DOUBLED. 
THE  USER  MAY  CHANGE  PRMT(4)  BY  MEANS  OF  HIS 
OUTPUT  SUBROUTINE. 

NO  INPUT  PARAMETER.  SUBROUTINE  DRKGS  INITIALIZES 
PRMT(5)=0.  IF  THE  USER  WANTS  TO  TERMINATE 
SUBROUTINE  DRKGS  AT  ANY  OUTPUT  POINT,  HE  HAS  TO 
CHANGE  PRMT (5)  TO  NON-ZERO  BY  MEANS  OF  SUBROUTINE 
OUTP.  FURTHER  COMPONENTS  OF  VECTOR  PRMT  ARE 
FEASIBLE  IF  ITS  DIMENSION  IS  DEFINED  GREATER 
THAN  5.  HOWEVER  SUBROUTINE  DRKGS  DOES  NOT  REQUIRE 
AND  CHANGE  THEM.  NEVERTHELESS  THEY  MAY  BE  USEFUL 
FOR  HANDING  RESULT  VALUES  TO  THE  MAIN  PROGRAM 
(CALLING  DRKGS)  WHICH  ARE  OBTAINED  BY  SPECIAL 
MANIPULATIONS  WITH  OUTPUT  DATA  IN  SUBROUTINE  OUTP. 
DOUBLE  PRECISION  INPUT  VECTOR  OF  INITIAL  VALUES 
(DESTROYED).  LATERON  Y  IS  THE  RESULTING  VECTOR  OF 
DEPENDENT  VARIABLES  COMPUTED  AT  INTERMEDIATE 
POINTS  X. 

DOUBLE  PRECISION  INPUT  VECTOR  OF  ERROR  WEIGHTS 
(DESTROYED).  THE  SUM  OF  ITS  COMPONENTS  MUST  BE 
EQUAL  TO  1.  LATERON  DERY  IS  THE  VECTOR  OF 
DERIVATIVES,  WHICH  BELONG  TO  FUNCTION  VALUES  Y  AT 
INTERMEDIATE  POINTS  X. 
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NDIM  -  AN  INPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 
EQUATIONS  IN  THE  SYSTEM. 

IHLF  -  AN  OUTPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 

BISECTIONS  OF  THE  INITIAL  INCREMENTe  IF  IHLF  GETS 
GREATER  THAN  10,  SUBROUTINE  DRKGS  RETURNS  WITH 
ERROR  MESSAGE  IHLF=11  INTO  MAIN  PROGRAM.  ERROR 
MESSAGE  I HL F= 1 2  OR  IHLF=13  APPEARS  IN  CASE 
PRMT(3)=0  OR  IN  CASE  SIGN (PRMT(3) ) .NE . SIGN (PRMT(2) - 
PRMT(l))  RESPECTIVELY. 

FCT  -  THE  NAME  OF  AN  EXTERNAL  SUBROUTINE  USED.  THIS 

SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDES  DERY  OF 
THE  SYSTEM  TO  GIVEN  VALUES  X  AND  Y.  ITS  PARAMETER 
LIST  MUST  BE  X,Y ,DERY .  SUBROUTINE  FCT  SHOULD 
NOT  DESTROY  X  AND  Y. 

OUTP  -  THE  NAME  OF  AN  EXTERNAL  OUTPUT  SUBROUTINE  USED. 

ITS  PARAMETER  LIST  MUST  BE  X,Y , DERY , IHLF, NDIM, PRMT. 
NONE  OF  THESE  PARAMETERS  (EXCEPT,  IF  NECESSARY, 
PRMT(4) ,PRMT (5) , . . . )  SHOULD  BE  CHANGED  BY 
SUBROUTINE  OUTP.  IF  PRMT(5)  IS  CHANGED  TO  NON-ZERO, 
SUBROUTINE  DRKGS  IS  TERMINATED. 

AUX  -  DOUBLE  PRECISION  AUXILIARY  STORAGE  ARRAY  WITH  8 
ROWS  AND  NDIM  COLUMNS. 

REMARKS 

THE  PROCEDURE  TERMINATES  AND  RETURNS  TO  CALLING  PROGRAM,  IF 

(1)  MORE  THAN  10  BISECTIONS  OF  THE  INITIAL  INCREMENT  ARE 
NECESSARY  TO  GET  SATISFACTORY  ACCURACY  (ERROR  MESSAGE 
I HL  F= 1 1 ) , 

(2)  INITIAL  INCREMENT  IS  EQUAL  TO  0  OR  HAS  WRONG  SIGN 
(ERROR  MESSAGES  IHLF=12  OR  IHLF=13), 

(3)  THE  WHOLE  INTEGRATION  INTERVAL  IS  WORKED  THROUGH, 

(4)  SUBROUTINE  OUTP  HAS  CHANGED  PRMT(5)  TO  NON-ZERO. 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
THE  EXTERNAL  SUBROUTINES  FCT(X,Y,DERY)  AND 
OUTP(X,Y, DERY, IHLF, NDIM, PRMT)  MUST  BE  FURNISHED  BY  THE  USER. 

METHOD 

EVALUATION  IS  DONE  BY  MEANS  OF  FOURTH  ORDER  RUNGE-KUTTA 
FORMULAE  IN  THE  MODIFICATION  DUE  TO  GILL.  ACCURACY  IS 
TESTED  COMPARING  THE  RESULTS  OF  THE  PROCEDURE  WITH  SINGLE 
AND  DOUBLE  INCREMENT. 

SUBROUTINE  DRKGS  AUTOMATICALLY  ADJUSTS  THE  INCREMENT  DURING 
THE  WHOLE  COMPUTATION  BY  HALVING  OR  DOUBLING.  IF  MORE  THAN 
10  BISECTIONS  OF  THE  INCREMENT  ARE  NECESSARY  TO  GET 
SATISFACTORY  ACCURACY,  THE  SUBROUTINE  RETURNS  WITH 
ERROR  MESSAGE  IHLF=11  INTO  MAIN  PROGRAM. 

TO  GET  FULL  FLEXIBILITY  IN  OUTPUT,  AN  OUTPUT  SUBROUTINE 
MUST  BE  FURNISHED  BY  THE  USER. 


SUBROUTINE  DRKGS(PRMT,Y ,DERY,NDIM,IHLF,FCT,OUTP,AUX) 


DIMENSION  Y(1 )  ,DERY(1 ) ,AUX(8,1 ) ,A(4) ,B(4) ,C(4) ,PRMT(1 ) 

DOUBLE  PRECISION  PRMT.Y.DERY ,AUX,A,B >C,X,XEND,H ,AJ ,BJ ,CJ ,R1  ,R2, 
1DELT 

DO  1  1  =  1  ,NDIM 

1  AUX(8,I )=.066666666666666667D0*DERY(I ) 

X=PRMT(1 ) 

XEND=PRMT(2) 

H=PRMT(3) 

PRMT(5)=0.D0 
CALL  FCT(X,Y,DERY) 

ERROR  TEST 

I F ( H* ( XEN  D-X ) ) 38 ,37,2 

PREPARATIONS  FOR  RUNGE-KUTTA  METHOD 

2  A(1 )=.5DO 

A(2)=. 29289321 881345248D0 
A  ( 3 )  =  1 .7071067811 865475D0 
A(4)  = .  1 6666666666666667D0 
B ( 1 )=2 . DO 
B { 2 ) = 1 .DO 
B(3)=l .DO 
B(4)=2 . DO 
C(1)=.5D0 

C(2)=. 29289321 881345248D0 
C ( 3 ) = 1 . 707106781 1 865475D0 
C(4)=.5D0 

PREPARATIONS  OF  FIRST  RUNGE-KUTTA  STEP 
DO  3.  1  =  1 ,  NDIM 
AUX ( 1 ,1)=Y(I) 

AUX(2,I)=DERY (I ) 

AUX  (3 ,1  )=0.D0 

3  AUX(6,I)=0.D0 
IREC=0 
H=H+H 
IHLF--1 
ISTEP=0 
IEND=0 


START  OF  A  RUNGE-KUTTA  STEP 

4  IF( (X+H-XEND)*H)7 ,6 ,5 

5  H=XEND-X 

6  IEND=1 


I  A.  '■  I,  S  3  T  r  1 1 


RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 

7  CALL  OUTP(X,Y,DERY,IREC,NDIM,PRMT) 

I F ( PRMT ( 5 ) )40 ,8 ,40 

8  ITEST=0 

9  I STEP= I STEP+ 1 


START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
J=1 

10  AJ=A( J) 

BJ=B(J) 

CJ=C(J) 

DO  11  1  =  1  ,NDIM 
R1=H*DERY(I ) 

R2=AJ*(R1 -BJ*AUX(6 ,1 ) ) 

Y  ( I  )=Y ( I )+R2 

R?=R?+R?+R? 

11  AUX(6 ,1 )=AUX(6 ,1 )+R2-CJ*Rl 
IF(J-4)12, 15,15 

12  J=J+1 
IF(J-3)13 ,14,13 

13  X=X+.5D0*H 

14  CALL  FCT(X,Y ,DERY ) 

GO  TO  10 

END  OF  INNERMOST  RUNGE-KUTTA  LOOP 


TEST  OF  ACCURACY 

15  IF(ITEST)! 6 ,1 6 ,20 

IN  CASE  ITEST=0  THERE  IS  NO  POSSIBILITY  FOR  TESTING  OF  ACCURACY 

16  DO  17  1=1,  NDIM 

17  AUX(4 ,1 )=Y ( I ) 

ITEST=1 

ISTEP=ISTEP+ISTEP-2 

18  IHLF=IHLF+1 
X=X-H 
H=.5D0*H 

DO  19  1=1 ,  NDIM 
Y ( I  )=AUX( 1 ,1 ) 

DERY ( I )=AUX(2 ,1 ) 

19  AUX(6,I)=AUX(3,I) 

GO  TO  9 

IN  CASE  I TEST= 1  TESTING  OF  ACCURACY  IS  POSSIBLE 

20  IMOD= I  STEP/2 

IF(ISTEP-I MOD- 1 MOD ) 2 1  ,23,21 
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21  CALL  FCT(X,Y,DERY) 

DO  22  1  =  1  ,NDIM 
AUX ( 5 ,1 )=Y (I ) 

22  AUX(7 ,1 )=DERY ( I ) 

GO  TO  9 

COMPUTATION  OF  TEST  VALUE  DELT 

23  DELT=OoDO 

DO  24  1  =  1  ,NDIM 

24  DELT=DELT+AUX(8,I)*DABS(AUX(4,I)-Y(I)) 
I F ( DELT-PRMT ( 4 ) ) 28 ,28 ,25 

ERROR  IS  TOO  GREAT 

25  I F( IHLF-1 0) 26 ,36,36 

26  DO  27  1  =  1  ,NDIM 

27  AUX ( 4,1 )=AUX(5 ,1 ) 

ISTEP=ISTEP+ISTEP-4 

X=X-H 

IEND=0 
GO  TO  18 

RESULT  VALUES  ARE  GOOD 

28  CALL  FCT(X,Y,DERY) 

DO  29  1  =  1  ,NDIM 
AUX ( 1 ,1 )=Y (I ) 

AUX( 2 ,1 )  =  DERY ( I ) 

AUX ( 3 ,1 )=AUX( 6 ,1 ) 

Y ( I )=AUX(5 ,1 ) 

29  DERY ( I )=AUX( 7 ,1 ) 

CALL  OUTP ( X-H  ,Y ,DERY ,IHLF,NDIM,PRMT ) 

I F ( P  RMT ( 5 ) ) 40, 30, 40 

30  DO  31  1  =  1  ,NDIM 
Y ( I )=AUX(1 ,1 ) 

31  DERY ( I )=AUX(2 ,1 ) 

I REC= I HL  F 

I F( IEND)32 ,32 ,39 

INCREMENT  GETS  DOUBLED 

32  I  HLF= I HL  F- 1 

I STEP= I  STEP/2 
H=H+H 

I F ( I HLF ) 4 ,33 ,33 

33  I MOD= I STEP/2 

IF( ISTEP-IM0D-IM0D)4 ,34,4 

34  IF(DELT- „02D0*PRMT(4) ) 35 ,35,4 

35  I HLF= IHLF-1 

I STEP= I  STEP/2 
H=H+H 
GO  TO  4 
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RETURNS  TO  CALLING  PROGRAM 

36  IHLF=1 1 

CALL  FCT(X,Y,DERY) 

GO  TO  39 

37  IHLF=1 2 
GO  TO  39 

38  IHLF=1 3 

39  CALL  OUTP(XsY  ,DERY ,IHLF,NDIM,PRMT) 

40  RETURN 
END 
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APPENDIX  B 


NUMERICAL  INTEGRATION 


B.l  Simpson's  Rule 

The  solution  of  the  system  of  differential  equations,  by  the 
scheme  described  in  Chapter  III,  gives  the  value  of  angle  $  at  equal 
intervals  of  s’.  From  this  equal-interval  data  "x  and  y  given  by  the 
expressions 


ft 

x  =  /  cos<f>  ds 
0 


ft 

-  _  /  sine})  ds 
y  ~  0 


B.l 


are  determined  by  the  application  of  Simpson's  Rule  (5). 

According  to  Simpson's  Rule,  the  vector  of  integral  values 


z.  =  z(xi)  =  /  y(x)  dx 


e  ®  o  3 


n) 


with 


x .  =  a  +  ( i  -  1 )  h 


is  given  by 


B.  2 
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To  use  this  formula,  minimum  three  points  are  required. 

B.2  Newton's  Three-Eighth  Rule 

The  major  disadvantage  of  Simpson's  Rule  is  that  is  requires 
the  function  values  at  odd  number  of  points.  This  is  overcome  by 
employing  it  in  combination  with  Newton's  Three-Eighth  Rule  (5),  which 
gives 


z*  =  z.0  +  -^r(y.o  +  3y.0  +  3y.  n+y.)  B.3 

J  J-3  8  ^  j-3  Jj-2  Jj-1  Jy 

B.3  Accuracy  of  the  Method 

All  the  calculations  are  carried  out  in  double-precision, 
that  is  carrying  16  significant  figures.  Therefore  it  is  reasonable 
to  assume  that  round  off  errors  are  negligible.  Truncation  error  is 

5 

of  the  order  h  in  both  methods.  The  value  of  h  used  in  the  present 
problem  is  of  the  order  10"^.  Therefore  truncation  error  too  is  negligible. 
A  Fortran  IV  Subroutine  DQSF,  based  on  the  method  described  above  is 
given  in  section  B.4. 
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B.4  SUBROUTINE  DQSF 


PURPOSE 

TO  COMPUTE  THE  VECTOR  OF  INTEGRAL  VALUES  FOR  A  GIVEN 
EQUIDISTANT  TABLE  OF  FUNCTION  VALUES. 


USAGE 

CALL  DQSF  (H,Y,Z,NDIM) 


DESCRIPTION  OF  PARAMETERS 

H  -  DOUBLE  PRECISION  INCREMENT  OF  ARGUMENT  VALUES. 

Y  -  DOUBLE  PRECISION  INPUT  VECTOR  OF  FUNCTION  VALUES. 

Z  -  RESULTING  DOUBLE  PRECISION  VECTOR  OF  INTEGRAL 

VALUES.  Z  MAY  BE  IDENTICAL  WITH  Y. 

NDIM  -  THE  DIMENSION  OF  VECTORS  Y  AND  Z. 


REMARKS 

NO  ACTION  IN  CASE  NDIM  LESS  THAN  3. 


SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


METHOD 

BEGINNING  WITH  Z(1)=0,  EVALUATION  OF  VECTOR  Z  IS  DONE  BY 
MEANS  OF  SIMPSONS  RULE  TOGETHER  WITH  NEWTONS  3/8  RULE  OR  A 
COMBINATION  OF  THESE  TWO  RULES.  TRUNCATION  ERROR  IS  OF 
ORDER  H**5  (I.E.  FOURTH  ORDER  METHOD). 

SUBROUTINE  DQSF(H,Y ,Z,NDIM) 


DIMENSION  Y( 1 ) ,  Z(l) 

DOUBLE  PRECISION  Y,Z,H,HT,  SUM!,  SUM2 ,  AUX,  AUX1 ,  AUX2 

HT= .33333333333333333D0*H 
I F ( NDIM- 5 ) 7,8,1 

NDIM  IS  GREATER  THAN  5.  PREPARATIONS  OF  INTEGRATION  LOOP 
1  SUM1 = Y (2)+Y(2) 

SUM1=SUM1+SUM1 
SUM1=HT*(Y ( 1 )+SUMl+Y(3) ) 

AUX1=Y (4)+Y (4) 

AUX1=AUX1+AUX1 

AUX1=SUM1+HT*(Y(3)+AUX1+Y (5) ) 

AUX2=HT*(Y(1 )+3. 875D0*(Y (2)+Y (5) )+2 „ 625D0*(Y(3)+Y(4) )+Y(6) ) 
SUM2=Y(5)+Y(5) 

SUM2=SUM2+SUM2 

SUM2=AUX2-HT*(Y(4)+SUM2+Y (6) ) 
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Z  ( 1 ) =0 . DO 
AUX=Y(3)+Y(3) 

AUX=AUX+AUX 

Z(2)=SUM2-HT*(Y (2)+AUX+Y(4) ) 

Z ( 3) =SUM1 
Z(4)=SUM2 
I F (NDIM-6) 5,5,2 
INTEGRATION  LOOP 

2  DO  4  I=7,NDIM,2 
S  UM1 =AUX 1 
SUM2=AUX2 

AUX1 =Y ( 1-1 )+Y(I-l ) 

AUX1=AUX1+AUX1 

AUX1 =SUM1 +HT*(Y( I -2)+AUXl+Y (I ) ) 

Z(I-2)=SUM1 
I F( I -NDI M) 3 ,6,6 

3  AUX2=Y(I)+Y(I) 

AUX2=AUX2+AUX2 

AUX2=SUM2+HT*(Y( I - 1 )+AUX2+Y ( 1+ 1 ) ) 

4  Z ( I  —  1 )=SUM2 

5  Z (NDI M- T )=AUX1 
Z(NDIM)=AUX2 
RETURN 

6  Z (NDI M- 1 )=SUM2 
Z(NDIH)=AUX1 
RETURN 

END  OF  INTEGRATION  LOOP 

7  I F  (NDI  M-3 )  12,11  ,8 

NDIM  IS  EQUAL  TO  4  OR  5 

8  SUM2=1.125D0*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4)) 
SUM1=Y(2)+Y(2) 

SUM1=SUM1+SUM1 
SUM1=HT*( Y( 1 )+SUMl+Y(3) ) 

Z(1)=0.D0 

AUX1=Y(3)+Y(3) 

AUX1=AUX1+AUX1 

Z(2)=SUH2-HT*(Y(2)+AUX1+Y  (4) ) 

IF(NDIM-5)10,9,9 

9  AUX1=Y (4)+Y(4) 

AUX1=AUX1+AUX1 

Z( 5)=SUM1+HT*(Y(3)+AUX1+Y(5) ) 

10  Z(3)=SUM1 

Z(4)=SUM2 
RETURN 
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NDIM  IS  EQUAL  TO  3 

11  SUM1 =HT*( 1 . 25D0*Y ( 1 )+Y ( 2 )+Y ( 2) - . 25D0*Y  ( 3) ) 
SUM2=Y(2)+Y(2) 

SUM2=SUM2+SUM2 
Z(3)=HT*(Y(1 )+SUM2+Y(3) ) 

Z ( 1 ) =0 . DO 
Z(2)=SUM1 

12  RETURN 
END 


